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Ab initio computational methods for electronic transport in nanoscaled systems are an invaluable 
tool for the design of quantum devices. We have developed a flexible and efficient algorithm for 
evaluating I-V characteristics of atomic junctions, which integrates the non- equilibrium Green's 
function method with density functional theory. This is currently implemented in the package 
Smeagol. The heart of Smeagol is our novel scheme for constructing the surface Green's functions 
describing the current/voltage probes. It consists of a direct summation of both open and closed 
scattering channels together with a regularization procedure of the Hamiltonian, and provides great 
improvements over standard recursive methods. In particular it allows us to tackle material systems 
with complicated electronic structures, such as magnetic transition metals. Here we present a 
detailed description of Smeagol together with an extensive range of applications relevant for the two 
burgeoning fields of spin and molecular-electronics. 



PACS numbers: 75.47. Jn, 72.10.Bg, 73.63.-b 



I. INTRODUCTION 

The study of electronic transport through devices com- 
prising only a handful of atoms is becoming one of the 
most fascinating branch of modern solid state physics. 
The field was initiated with the advent of the scanning 
tunnelling microscope (STM)i and at present comprises 
a multitude of applications which span over several disci- 
plines and encompass different technologies, from build- 
ing blocks for revolutionary computer architectures, to 
disposable electronics, to diagnostic tools for genetically 
driven medicine. Clearly many of these devices will soon 
change and enhance the quality of our daily life. 

Fully functional molecular memories 2 and logic 
gates 3 ^ have been already demonstrated suggesting a 
possible roadmap to the post-silicon era. These should 
produce future generations of computers, together with 
magnetic data storage devices exceeding the Terabit /in 2 
storage limit. The readout of such high-density data stor- 
age media will be achieved using nanoscale devices with 
magnetic atomic point contacts^. 

At the same time hybrid molecular devices are becom- 
ing increasingly popular in multifunctional sensor design, 
demonstrating a sensitivity orders of magnitude supe- 
rior to that achievable with conventional methods. These 
molecular devices include for example carbon nanotubes 
detectors for NO^ and nerve agents 8 -, nanowire-based 
virus detectors 9 - and chemical sensors^. The near future 
should see the development of on-chip nanolabs able to 
sense a particular signature of gene or protein expres- 
sion and therefore be able to diagnose various diseases. 
These will be formidable tools for the study of biological 
systems and in the field of preventive medicineii. 

In addition to this large experimental activity an 
equally large effort has been devoted to the develop- 
ment of efficient computational methods for evaluating 
I-V characteristics of nanoscale devices. This is quite 



a remarkable theoretical challenge since advanced quan- 
tum transport algorithms must be combined with state of 
the art electronic structure methods. Ideally these tools 
should be able to include strong correlation as well as 
inelastic effects, and they should be suitable for describ- 
ing large systems (easily scalable methods). Furthermore 
in order to compare directly to experiments the detailed 
knowledge of the atomic configuration is needed. 

Modern theory of quantum transport has developed a 
range of methods for calculating transport in nanoscale 
conductors. Broadly speaking these can be divided into 
two main classes: 1) steady state algorithms, and 2) 
time dependent schemes. The first are based upon the 
assumption that, regardless of the details of a possible 
transient, a steady state is always achieved. The current 
through the entire device is calculated as a balance of 
currents entering and leaving a given scattering region, 
either using scattering theorjii2ii 3 rfi2ii£*i 6 ^i 7 *i&, or by solv- 
ing a master equation 19,20,21 . A multitude of variations 
over this generic scheme are available 2 ^, depending on 
the underlying assumption leading to the steady state, 
the details of the electronic structure method employed, 
and the way in which the external potential is introduced 
in the calculation. Interestingly most of the methods 
can be demonstrated to be applicable to cases of non- 
interacting electrons 2 ^, although their equivalence is not 
demonstrated for the interacting case. Among these algo- 
rithms a particular place is occupied by implementations 
of the non-equilibrium Green function (NEGF) 12 : 13 : 14 : 15 
method within density functional theory (DFT)—^. 
This approach, which is based on equilibrium DFT to 
describe the electronic structures, has the advantage of 
being conceptually simple, and computationally easy and 
versatile to implement 26 ' 27 : 28 ' 29 : 30 . 

Time dependent methods are at an earlier stage of de- 
velopment. These investigate the time-evolution of the 
electronic charge density of a system, brought out of equi- 
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librium by a time-dependent perturbing potential. To 
the best of our knowledge two fundamentally different 
methods have been proposed to date. The first considers 
infinite non-periodic systems, with an external potential 
introduced as a time dependent perturbationSi^i. The 
time-evolution of the density matrix is studied with time- 
dependent density functional theory (TDDFT)? 2 i??. An 
alternative approach consists of placing the system of in- 
terest in a large capacitor. Such a capacitor is charged 
at t=0 and the time-dependent de-charging process is 
investigated^*^. Generally speaking these methods are 
computationally intensive, since the need to perform the 
time evolution adds to the computational overheads of 
standard static schemes. However they should be able 
to address transport limits such as Coulomb blockade 
or resonant tunnelling, which are otherwise difficult to 
describe. Interestingly, important information can be 
extracted from the static limit of the time-dependent 
problemSi*i2£, and this can help in designing more ac- 
curate static methods and in understanding their limita- 
tions. 



II. NON-EQUILIBRIUM TRANSPORT 
METHOD 

In this section we describe in details our compu- 
tational technique. The underlying assumption used 
throughout this work is that all the quantities associ- 
ated to the electronic structure (Hamiltonian, density 
matrix, Green's functions, etc.) can be written over a 
localized atomic orbital (LAO) basis set of some kind 
ipin = ip^(f — Ri), where Ri is the position of the i- 
th nucleus and fi — n,l,m is a collective index span- 
ning, the angular momentum (I, m) and the orbital n. 
Note that in general the index n can run over different 
radial functions corresponding to the same angular mo- 
mentum, according to the multiple-zetas schemed. In 
this way a generic operator O is represented by a finite 
N x N matrix (N is the total number of degrees of free- 
dom in the system) whose matrix elements are simply 
Oi fj,, jv Note also that the functions tpi^ are generally 
non-orthogonal and the overlap matrix S is defined as 
Si* i» = 1 1>%?- Ri)Mr~ Rj) d 3 f. 



Here we present in details our recently developed quan- 
tum transport code SmeagoS^. Smeagol is a DFT im- 
plementation of the non-equilibrium Green's function 
method, which has been specifically designed for mag- 
netic materials. The main core of Smeagol is our original 
technique for constructing the leads self-energies^ which 
avoids the standard problems of recursive methods 12 
and allows us to describe devices having current /voltage 
probes with a complicated electronic structure. In ad- 
dition Smeagol has been constructed to be a modular 
and scalable code, with particular emphasis on heavy 
parallelization, to facilitate large scale simulations. In 
its present form Smeagol is parallel over fc-space, real 
space and energy and furthermore it can deal with spin- 
polarised systems, including spin non-collinearity. A par- 
tial description of the code has already been provided 4 ^, 
which should be incorporated with this more detailed de- 
scription. 



The paper is organized as follows. In the next section 
we introduce our method and its main technical imple- 
mentations. In particular we set the problem, explain 
how to construct the leads self-energies, and describe the 
strategy used for calculating the electrostatic potential. 
Then we present a series of calculations for systems rel- 
evant to either spin- or molecular-electronics. These ad- 
dress specific aspects of Smeagol such as the electrostat- 
ics, the spin polarization and the spin non-collinearity. 
Finally in the appendices we describe in more details the 
self-energy algorithm, we recall the theoretical founda- 
tions of the NEGF formalism, and we establish a connec- 
tion with TDDFT. 



A. Problem setup 

Smeagol has been designed to describe two-terminal 
conductance experiments, where two current /voltage 
electrodes of macroscopic size sandwich a nanometer- 
sized device (a molecule, an atomic point contact, a tun- 
neling barrier...). Let us present the problem from three 
different perspectives: the thermodynamics, the Hamil- 
tonian and the electrostatics (see figure 

From a thermodynamic point of view the system is 
modelled two bulk leads and a central region. The lat- 
ter includes the actual device and, for reasons that will 
be clear later, part of the leads. Therefore we call such 
a central region an "extended molecule" (EM). The two 
current/ voltage leads are kept at two different chemical 
potentials respectively /il and /xr and are able to ex- 
change electrons with the EM. Note that when the ap- 
plied bias is zero (/xl = (J-r), this system of interacting 
electrons is in thermodynamic equilibrium and may be 
regarded as a grand canonical ensemble. When the bias 
is applied however /il 7^ A*r and the current will flow. 
Then the prescription for establishing the steady state is 
that of adiabatically switching on the coupling between 
the leads and the EM 14 i 42 i 43 . 

At the Hamiltonian level the system under investiga- 
tion is described by an infinite hermitian matrix TL. This 
however has a rather regular structure. First notice that 
the two semi-infinite current /voltage probes are defect- 
free crystalline metals. These have a regular periodic 
structure and a unit cell along which the direction of the 
transport can be defined. At this point it is important 
to notice that because of the LAO basis set this matrix 
will be rather sparse. It is convenient to introduce the 
concept of principal layer (PL). A principal layer is the 
smallest cell that repeats periodically in the direction of 
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Figure 1: Schematic two terminal device, (a) Thermodynam- 
ical aspect: two leads are kept respectively at the chemical 
potentials ^ and and are able to exchange electrons with 
the central region (extended molecule EM), (b) Hamiltonian 
description: two block diagonal infinite matrices describe the 
semi-infinite current /voltage probes, and a finite matrix Hm 
describes the extended molecule. Hq is a finite Hamiltonian 
matrix describing one principal layer, while Hi describes the 
interaction between two adjacent principal layers, (c) Electro- 
statics: the two current /volt age leads have a constant average 
potential /^l/r = Ep ± V/2, and the potential drop occurs 
within the extended molecule. 



the transport constructed in such a way to interact only 
with the nearest neighbour PLs. This means that all the 
matrix elements between atoms belonging to two non- 
adjacent PLs vanish. For example take a linear chain of 
hydrogen atoms described by a nearest neighbour tight- 
binding model then one atom forms the unit cell. How- 
ever if nearest and next nearest neighbour elements are 
included then the PL will contain two atoms etc (for ex- 
amples see appendix 0. 

We then define Hq as the N x N matrix describing all 
interactions within a PL, where N is the total number 
of degrees of freedom (basis functions) in the PL (note 
that we use calligraphic symbols - TL - for infinitely- 
dimensional matrices and capitalized letters - H - for 
finite matrices). Similarly Hi is the N x N matrix de- 
scribing the interaction between two PLs. Finally Hm is 
the M x M matrix describing the extended molecule and 
Hlm (-Srm) is the N x M matrix containing the interac- 
tion between the last PL of the left-hand side (right-hand 
side) lead and the extended molecule. The final form of 
TL is 
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H-i = Hi, H M l = h[ m and H M r = #r M - In tnis form 
TL has the same structure as the Hamiltonian of a one- 
dimensional system as shown in figure^}. However this 
is not the most general situation and does not apply if a 
magnetic field is present for example. 

Note that the overlap matrix S has exactly the same 
structure of TL. Therefore we adopt the notation So, Si, 
Slm, Srm and Sm for the various blocks of S, in com- 
plete analogy with their Hamiltonian counterparts. Here 
the principal layer, defined by TL is used for both the iS 
and the TL matrix, even though the range of S can be 
considerably shorter than that of TL. 

Let us now discuss the electrostatics of the problem 
(figure ^;). The main consideration here is that the 
current/ voltage probes are made from good metals and 
therefore preserve local charge neutrality. For this rea- 
son the effect of an external bias voltage on the leads will 
produce a rigid shift of the whole spectrum, i.e. of all the 
on-site energies. In contrast a non-trivial potential profile 
will develop over the extended molecule, which needs to 
be calculated self-consistently. Importantly the resulting 
self-consistent electrostatic potential must match that of 
the leads at the boundaries of the EM. If this does not 
happen, the potential profile will develop a discontinuity 
with the generation of spurious scattering. Therefore, in 
order to achieve a good match of the electrostatic po- 
tential, several layers of the leads are usually included 
in the extended molecule. Their number ultimately de- 
pends upon the screening length of the leads, but in most 
situations a few (between two and four) atomic planes are 
sufficient. 

Even in the case of extremely short screening length, 
it is however good practice to include a few planes of the 
leads in the extended molecule because the electrodes 
generally have reconstructed surfaces, which might un- 
dergo additional geometrical reconstructions when bond- 
ing to the nanoscale device (e.g. molecules attached to 
metallic surfaces through corrosive chemical groups). 

We conclude this section with some comments about 
the application of periodic boundary conditions in the 
transverse direction perpendicular to that of the trans- 
port. The setup of a typical experiment is that of two 
very large current voltage probes sandwiching a tiny re- 
gion which is responsible for most of the resistance. The 
ideal description would be that of two infinite leads (with 
infinite cross sections) and a finite scattering region. Un- 
fortunately this problem is intractable since both Ho and 
Hi become infinitely-dimensional. Therefore one has to 
consider some approximations. 

The first option is to use leads with finite cross section. 
In this case, no periodic boundary conditions are required 
and the whole system is quasi-one-dimensional. However 
special care must be taken when choosing the cross sec- 
tion of the leads in order to avoid quantum confinement 
effects. It is also worth noting that leads with very small 
cross section make the use of the Landauer formula for 
transport^ questionable. As a rule of thumb the linear 
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dimension of the cross section should be several times the 
Fermi wavelength of the material forming the leads and 
there should be several open scattering channels. 

The second option is to use periodic boundary condi- 
tions. In this case the system is repeated periodically 
in the transverse direction, meaning that the extended 
molecule is also repeated periodically. Clearly quan- 
tum confinement effects are eliminated, but one should 
be particularly careful in order to eliminate the spuri- 
ous interaction between the mirror images of the ex- 
tended molecule. Therefore large unit cells must be 
employed even when periodic boundary conditions are 
used. However from a formal point of view the use of 
periodic boundary conditions does not change the prob- 
lem setup. All the matrices (Ho, Hi etc.) now depend 
on the transverse k- vector used, and the infinite prob- 
lem transforms in a collection of /c-dependent quasi one- 
dimensional problems. This dependence is implicitly as- 
sumed whenever necessary throughout the rest of the pa- 
per. 



B. Green's functions for an open system 

We are dealing with an infinite-dimensional hermitian 
problem, which is intractable, because the wave-functions 
deep inside the leads have a plane-wave form. These 
can be calculated by computing the band-structure of 
an infinite chain of PL's. The main computational ef- 
fort is therefore focussed upon the problem of describing 
the scattering of plane- waves from one lead to the other 
across the EM. The problem is solved by computing the 
retarded Green's function C/ R for the whole system by 
solving the Green's function equation 



[e+S-H] G R (E) = 2, 



(2) 



where I is an infinitely-dimensional identity matrix, e + = 
lim,5^o+ E + iS and E is the energy. The same equation 
explicitly using the block-diagonal structure of both the 
Hamiltonian and the overlap matrix (we drop the symbol 
"R" indicating the retarded quantities) is of the form 



J 



— Hl £ + <Slm — H.LM 

e + <?ML — Hml £ + Sm — Hm £ + Smr. — Hmr 
e + 6>R,M — 7~Lrm £ + 5r — Hr, 



Gl 


Glm 


Glr. i 









Gml 


Gm 


Gmr. 


-: 




s 


Grl 


Grm 


Gr J 










(3) 



r 



where we have partitioned the Green's functions G into 
the infinite blocks describing the left- and right-hand side 
leads Gl and Gr, those describing the interaction between 
the leads and extended molecule Glm, Grm, the direct 
scattering between the leads Glr, and the finite block 
describing the extended molecule Gm- We have also in- 
troduced the matrices Hl, Hr, Hlm, Hrm and their 
corresponding overlap matrix blocks, indicating respec- 
tively the left- and right-hand side leads Hamiltonian and 
the coupling matrix between the leads and the extended 
molecule. Hm is an M x M matrix and Im is the M x M 
unit matrix. The infinite matrices, Hl and Hr describe 
the leads and have the following block-diagonal form 



H L = 
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(4) 



with similar expressions for Hr and the overlap S matrix 
counterparts. In contrast the coupling matrices between 
the leads and the extended molecule are infinite dimen- 
sional matrices whose elements are all zero except for a 
rectangular block coupling the last PL of the leads and 
the extended molecule. For example we have 



H 



LM 




(5) 



The crucial step in solving equation Jljjl is to write 
down the corresponding equation for the Green's func- 
tion involving the EM and surface PL's of the left and 
right leads and then evaluate the retarded Green's func- 
tion for the extended molecule G M . This has the form 



G M (E) = [e+S M -H M - T*(E) - " 



(6) 



where we have introduced the retarded self-energies for 
the left- and right-hand side lead 

E R (£) = (e+Sml - H ML )G L R (E) (e+S LM - H LM ) (7) 

and 

^r(^) — ( e+ 5 , MR — Hmr)Gr 1 (E) (£ + 5rm — Hrm) ■ (8) 

Here G L R and G™ are the retarded surface Green's 
function of the leads, i.e. the leads retarded Green's 
functions evaluated at the PL neighboring the extended 
molecule. Formally G^ (G R R ) corresponds to the right 
lower (left higher) block of the retarded Green's function 
for the whole left-hand side (right-hand side) lead. These 
are simply 



and 



G L ,n (E) = [e+S L -HLY 



G° R R (E) = [e+S R -HRY 



(9) 



(10) 
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Note that £° R (£° R ) is not the same as £? R defined 
in equation (J3J). In fact the former are the Green's func- 
tions for the semi-infinite leads in isolation, while the 
latter are the same quantities for the leads attached to 
the scattering region. Importantly one does not need to 
solve the equations @ and (|10(l for calculating the loads 
surface Green's functions and a closed form avoiding the 
inversion of infinite matrices can be provideo— . This will 
be discussed in what follows and in appendix 1X1 

Let us conclude this section with a few comments on 
the results obtained. The retarded Green's function G^ 
contains all the information about the electronic struc- 
ture of the extended molecule attached to the leads. In 
its close form given by the equation J§} it is simply the re- 
tarded Green's function associated to the effective Hamil- 
tonian matrix H e g 

H cS = H M + + £*(£) . (11) 

Note that H e g is not hermitian since the self-energies are 
not hermitian matrices. This means the the number of 
particles in the extended molecule is not conserved, as ex- 
pected by the presence of the leads. Moreover, since G^ 
contains all the information about the electronic struc- 
ture of the extended molecule in equilibrium with the 
leads, it can be directly used for extracting the zero-bias 
conductance G of the system. In fact one can simply 
apply the Fisher-Lee— 4& relation and obtain 

G = ^Tr[r L G$ r R G R ] , (12) 

where 

r a (£)=iK(£)-E;(E)t]. (13) 

In equation l|12f) all the quantities are evaluated at the 
Fermi energy E F . Clearly Tr[r L G$ r R G^] (E) is simply 
the energy dependent total transmission coefficient T(E) 
of standard scattering theory— 

Finally note that what we have elaborated so far is 
an alternative way of solving a scattering problem. In 
standard scattering theory one first computes the asymp- 
totic current carrying states deep into the leads (scatter- 
ing channels) and then evaluate the quantum mechanical 
probabilities for these channels to be reflected and trans- 
mitted through the extended molecule— In this case 
the details of the scattering region are often reduced to 
a matrix describing the effective coupling between the 
two surface PLs of the leads— In contrast the use of 
(|12|l describes an alternative though equivalent approach, 
in which the leads are projected out to yield a reduced 
matrix describing an effective EM. The current through 
surface PL's perpendicular to the transport direction are 
the same^i, the two approaches are equivalent and there 
is no clear advantage in using either one or the other. 
However, when the Hamiltonian matrix of the scatter- 
ing region Hyi is not known a priori, then the NEGF 
method offers a simple way of setting up a self-consistent 
procedure. 



C. Steady-state and self-consistent procedure 

Consider now the case in which the matrix elements of 
the Hamiltonian of the system are not known explicitly, 
but only their functional dependence upon the charge 
density p, H = T~L[p], is known. This is the most common 
case in standard mean field electronic structure theory, 
such as DFT. If no external bias is applied to the device 
(linear response limit) the Hamiltonian of the system can 
be simply obtained from a standard equilibrium DFT cal- 
culation and the procedure described in the previous sec- 
tion can be applied without any modification. However, 
when an external bias V is applied, the charge distri- 
bution of the extended molecule will differ from that at 
equilibrium since both the net charge and the electrical 
polarization are affected by the bias. This will deter- 
mine a new electrostatic potential profile with different 
scattering properties. 

These modifications will affect only the extended 
molecule, since our leads are required to preserve local 
charge neutrality. This means that the charge density 
and therefore the Hamiltonian of the leads are not mod- 
ified by the external bias applied. As discussed at the 
beginning the only effect of the external bias over the 
current/ voltage electrodes is that of a rigid shift of the 
on-site energies. The Hamiltonian then takes the form 

/H L +leV/2 H LM \ 

T~i = Hml Hm Hmr , (14) 

V Hbm H R -leV/2 J 

Note that the coupling matrices between the leads and 
the extended molecule are also not modified by the exter- 
nal bias, since by construction the charge density in the 
surface planes of the extended molecule matches exactly 
that of the leads. 

The Hamiltonian of the extended molecule 

H M =H M [p] (15) 

depends on the density matrix, which is calculated using 
the lesser Green's function G^ 12 i 13 i 14 i 15 i 42 i 43 

P=^-JdEG<(E), (16) 

so a procedure must be devised to compute this quantity. 

In equilibrium, G<(E) = -2ilm [G R (E)] f(E - p), 
so it is only necessary to consider the retarded Green's 
function, given by equation Alternatively, G R may 
be obtained from the eigenvectors of TL. 

Out of equilibrium, however, the presence of the 
leads establishes a non-equilibrium population in the 
extended molecule and G < is no longer equal to 
— 2ilm [G R ] f(E — p). The non-equilibrium Green's 
function formalism— iiiiiiS^SiMi provides the correct ex- 
pression (see appendix IBl: 

G< (E) = iG*(E)[r h f(E - p L ) + T R f(E - p R )]G^(E) 

(17) 
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where Ml/h, = cV/2, /(a;) is the Fermi function for a 
given temperature T, 

r L/R = r L/R (£±eV72) (is) 

and G^(-E) is given again by the equation © where now 
we replace £ L/R (-E) with £ L/R = £ l /r(£ ± eV/2). 

Finally the self consistent procedure is as follows. First 
a trial charge density p° is used to compute Hm from 
equation l(T5|l . Then Ll, L r and G^ are calculated from 
equations i|18fl . and These quantities are used to 

compute the Gj^ of equation (fT7|l . which is fed back in 
equation (|16|l to find a new density p . This process is 
iterated until a self-consistent solution is achieved, which 
is when \\p> — <C 1. At this point the problem is 

identical to that solved in the previous section (since the 
whole TL is now determined) and the current / can be 
calculated using24 

I=lJdE Tr[r L G^r R G^] [/ (E - ML ) - / (E - m )} . 

(19) 

Note that now the transmission coefficient depends on 
both the energy E and the bias V. 

Let us conclude this section with a note on how to 
perform the integrals of equations l|16l) and (|19l) . The one 
for the current is trivial since the two Fermi functions 
effectively cut the integration to give a narrow energy 
window between the chemical potentials of the leads. In 
addition the transmission coefficient, with the exception 
of some tunneling situations, is usually a smooth function 
of the energy. 

In contrast the integration leading to the density ma- 
trix (|16f) is more difficult, since the integral is unbound 
and the Green's function has poles over the real en- 
ergy axis. This however can be drastically simplified by 
adding and subtracting the term G^T R G^/ [E — p^) 
and by re-writing the integral Q16[) as the sum of two 
contributions p = p cq + pv 

Pcq = -^ J dEIm[G*]f(E-p L ), (20) 

and 

pv = hj dE G M r R G M if ( E - Mr) -f(E- Ml)] . 

(21) 

p eq can be interpreted as the density matrix at equilib- 
rium, i.e. the one obtained when both the reservoirs have 
the same chemical potential /iL, while pv contains all the 
corrections due to the non-equilibrium conditions. Com- 
putationally pv is bound by the two Fermi functions of 
the leads, as for the current J, and therefore one needs to 
perform the integration only in the energy range between 
the two chemical potentials. In contrast p oq is unbound, 
but the integral can be performed in the complex plane 
using a standard contour integral technique^, since G^ 
is both analytical and smooth. 



D. Surface Green's functions 

Let us now return to the question of how to calculate 
the self-energies for the leads. From the equations (7J and 
© it is clear that the problem is reduced to that of com- 
puting the retarded surface Green functions for the left- 
(G^ 11 ) and right-hand side (G R R ) lead respectively. This 
does not require any self-consistent procedure since the 
Hamiltonian is known and it is equal to that of the bulk 
leads plus a rigid shift of the on-site energies. However 
the calculation should be repeated several times since the 
E's depend both on the energy and the /c-vector. There- 
fore it is crucial to have a stable algorithm. 

There are a number of techniques in the literature to 
calculate the surface Green's functions of a semi-infinite 
system. These range from recursive methodsi2^ to semi- 
analytical constructions' 3 -. In Smeagol we have general- 
ized the scheme introduced by Sanvito et. al. 39 to non- 
orthogonal basis sets. This method gives us a prescrip- 
tion for calculating the retarded surface Green's function 
exactly. The main idea is to construct the Green's func- 
tion for an infinite system as a summation of Bloch states 
with both real and imaginary wave-vectors, and then to 
apply the appropriate boundary conditions to obtain the 
Green's function for a semi-infinite lead. 
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Figure 2: Infinite periodic system used as current/voltage 
probe and schematic diagram of the Hamiltonian. Ho and So 
are the matrices describing the Hamiltonian and the overlap 
within a PL, while Hi and Si are the same quantities cal- 
culated between two adjacent PLs. The arrow indicates the 
direction of transport. 

As explained above the Hamiltonian and the overlap 
matrices are arranged in a tridiagonal block form, having 
respectively Ho and So on the diagonal, and H\ and S± as 
the first off diagonal blocks (see figure |2J • Since we are 
dealing with an infinite periodic quasi one-dimensional 
system, the Schrodinger equation can be solved for Bloch 
states 

iP z = n^V^ (22) 

and reads 

[K + K ie lk + K_ ie - lk ] (j> k = , (23) 

where z — aoj with j integer and ao the separation be- 
tween principal layers, k is the wave- vector along the 
direction of transport (in units of n/ao), 4>k is a N- 
dimensional column vector and rik a normalization factor. 
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Here we introduce the N x N matrices 

Kq = H — ESo , 
K X =H X - ES l , 
K-y = B-x - ES-i . 



(24) 
(25) 
(26) 



Since the Green's functions are constructed at a given 
energy our task is to compute k(E) (both real and com- 
plex) instead of E(k) as conventionally done in band the- 
ory. A numerically efficient method to solve the "inverse" 
secular equation k = k(E) is to map it onto an equivalent 
eigenvalue problem. It is simple to demonstrate^ that 
the eigenvalues of the following 2N x 2N matrix 







(27) 



are e %k and that the upper N component of the eigenvec- 
tors are the vectors 4>k- Clearly for the solution of this 
eigenvalue problem one needs to invert K\. However, 
since K i is determined by the details of the physical sys- 
tem, the choice of basis set and of principal layer may be 
singular or severely ill-conditioned. This problem often 
originates from the fact that a few states within a PL do 
not couple to states in the nearest-neighbouring PLs, but 
it can also be due to the symmetry of the problem. For 
example in the case of ab initio derived matrices this be- 
comes unavoidable when one considers transition metals, 
where the strongly localized d shells coexist with rather 
delocalised s electrons. A possible solution to this prob- 
lem is to consider an equivalent generalized eigenvalue 
problem, which does not require matrix inversion. How- 
ever this solution is not satisfactory for two reasons. First 
the matrices still remain ill-conditioned and the general 
algorithm is rather unstable. Secondly for extreme cases 
we have discovered that the generalized eigenvalue solver 
cannot return meaningful eigenvalues (divisions by zero 
are encountered). We therefore decide to use an alter- 
native approach constructing a regularization procedure 
for eliminating the singularities of K\. This must be 
performed before starting the actual calculation of the 
Green's functions. We will return on this aspect in ap- 
pendix^ For the moment we assume that K\ has been 
regularized and it is neither singular nor ill-conditioned. 

When using orthogonal basis sets the knowledge of k 
and {<f>k} is sufficient to construct the retarded Green's 
function for the doubly-infinite system, which has the 
forrr*22i 



G zz , = 



Hz-^slv- 1 z>z' 



Yf^eM'-'^&V- 1 z<z> ' 



(28) 



where the summation runs over both real and imaginary 
ki. In equation l|28l) ki (fcj) are chosen to be the right- 
moving or right-decaying (left-moving or left-decaying) 
Bloch states, i.e. those with either positive group velocity 
or having /c-vector with positive imaginary part (negative 
group velocity or negative imaginary part), {(j)^ } are the 



corresponding vectors, and V is defined in reference |3£ 
Finally is just the dual of obtained from 



'LP*- 



(29) 
(30) 
(31) 



In the case of a non-orthogonal basis set the same ex- 
pression is still valid if V is now defined as follows 



N 



V 



(Hi ESl 



-iki it 



_ . (32) 

Finally the surface Green's functions for a semi-infinite 
system can be obtained from those of the doubly-infinite 
one by an appropriate choice of boundary conditions. For 
instance if we subtract the term 

N 

A z (z' z ) =^^ h e^^-^4 h ^e*«( z °- z ')4 i y-i , 

Lh 

(33) 

from G zz i of equation l|28l) we obtain a new retarded 
Green's function vanishing at z = zq. Note that 
A z (z! — zq) is a linear combination of eigenvectors and 
therefore does not alter the causality of G. 

In this way we obtain the final expression for the re- 
tarded surface Green's functions of both the left- and 
right-hand side lead 



r,0R _ 
u r — 



Lh 



Lh 



y-\(34) 



y- x .(35) 



These need to be computed at the beginning of the cal- 
culation only. 



E. DFT implementation and electrostatics 

The formalism presented in sections 111 Al through 111 Dl 
is rather general and is not specific of a particular func- 
tional dependence of the Hamiltonian upon the charge 
density. Therefore one can use on the same footing 
Hamiltonian theories ranging from parameterized self- 
consistent tight-binding methods^ to density functional 
theory 24-25 . Smeagol uses DFT as its main electronic 
structure method. 

At this point it is important to observe that Smeagol, 
and indeed any other NEGF DFT-based scheme, sim- 
ply uses the Kohn-Sham Hamiltonian 2 ^ as a single parti- 
cle Hamiltonian. This means that the non-equilibrium 
charge density obtained through the NEGF method 
(equations (|20|) and $21$ ) is not by any mean associ- 
ated with any variational principle and certainly does not 
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minimise the density functional, nor makes it stationary. 
The only exception is for zero bias, where the method 
presented here is just a clever alternative for solving an 
equilibrium problem for an infinite non-periodic system. 
Although it is common practice, it is therefore mislead- 
ing and incorrect to refer to our method as DFT-based 
NEGF, since the Hohenberg-Kohn theorem cannot be 
applied 24 . Some additional discussion over this issue can 
be found in reference [5l| . 

Although Smeagol is constructed in a simple and mod- 
ular way and can be readily interfaced with any DFT 
package based on a LAO basis set, for the present im- 
plementation we have used the existing code Siesta^. 
Siesta is a mature numerical implementation of DFT, 
which has been specifically designed for tackling prob- 
lems involving a large number of atoms. It uses norm 
conserving pseudopotentials in the separate Kleinman- 
Bylander form—, and most importantly a very efficient 
LAO basis set 41 i 54 i 55 . 

One important aspect that deserves a mention is the 
way in which we calculate the Hartree (electrostatic) po- 
tential for the extended molecule under bias. Clearly the 
easier and more transparent way would be that of solv- 
ing the Poisson's equation in real space with appropriate 
boundary conditions. However this usually is numeri- 
cally less efficient then solving it in fc-space using the 
fast Fourier transform algorithm. Siesta uses this second 
strategy and so does Smeagol. In Smeagol the electro- 
static potential is then calculated for the infinite system 
obtained by repeating periodically the extended molecule 
along the transport direction (see figure[3J). However, be- 
fore solving the Poisson's equation for such a system we 
add to the Hartree potential a saw-like term, whose drop 
is identical to the bias applied. For convergence reasons 
we often add at both edges of the scattering region two 
buffer layers, in which the external potential is only a 
constant and the density matrix is that of the leads and 
is not evaluated self-consistently. 

In summary a typical Smeagol calculation proceeds as 
follows. First one computes the leads self-energies over 
a range of energies E as large as the bandwidth of the 
materials forming the leads. These are then stored either 
in memory or on disk depending on the size of the leads. 
Then the proper Smeagol calculation is performed follow- 
ing the prescription described in the previous sections. 



III. TEST CASES 

We now present several test cases demonstrating the 
capabilities of Smeagol. They address key aspects of 
the code such as the electrostatics, the calculation of 
the transmission coefficients, the calculation of the I- 
V characteristic, the spin-polarization and the spin non- 
collinearity. 




Figure 3: Schematic representation of the electrostatic prob- 
lem. The real system (a) of an extended molecule sandwiched 
between two leads is mapped onto a fictitious periodic system 
(b), obtained by repeating the extended molecule in the direc- 
tion of the transport. The crucial point is that the potential 
profile in the unit cell of the periodic system is identical to 
that of the actual structure. 



A. Electrostatics: parallel plate capacitor 

As first simple test we present the case of a par- 
allel plate capacitor, constructed from two infinite fee 
gold surfaces separated by a vacuum region 12.3 A long. 
Clearly we do not expect transport across this device 
(with the exception of a tiny tunneling current), but it is 
a good benchmark of the Smeagol ability to describe the 
electrostatics of a device. 

The two gold surfaces are oriented along the (100) di- 
rection and the unit cell has only one atom in the cross 
section. The extended molecule comprises seven atomic 
planes in the direction of the transport, which is enough 
for achieving a good convergence of the Hartree potential 
(the Thomas-Fermi screening length in gold is ~0.6AS). 
For the calculation we use 100 fc-points in the full Bril- 
louin zone in the transverse direction, a single zeta basis 
set for the s, p and d orbitals and standard local density 
approximation (LDA) of the exchange and correlation 
potential. 

In figure^we present the planar average of the Hartree 
(electrostatic) potential Vh, the difference between the 
planar average of Hartree potential at finite bias and that 
at zero bias AV, and the difference Ap between of the 
planar average of the charge density along the direction of 
the transport for a given bias and that at zero bias. The 
quantities shown in the picture are those expected from 
the classical physics of a parallel plate capacitor. In the 
leads the electrostatic potential shows oscillations with a 
period corresponding to that of the separation between 
the gold planes, but with a constant average. In con- 
trast in the vacuum region the potential is much higher, 



9 




3 1 
o 

> n 



< 



1 


1 Volt 

2 Volts 










1 




i 




Figure 4: a) Planar average of the Hartree potential Vh for 
an infinite parallel-plate capacitor, b) Difference between the 
planar average of the Hartree potential at a given bias and 
that at zero bias AV. c) Difference Ap between of the planar 
average of the charge density along the direction of the trans- 
port for a given bias and that at zero bias. The dots indicates 
the position of gold planes. 



since there are no contributions from the nuclei, but it is 
uniform. If we eliminate the oscillations, by subtracting 
the zero bias potential from that obtained at finite bias 
(figure Et)) we obtain a constant potential profile in the 
leads and a linear drop in the vacuum region. Finally the 
macroscopic average of the charge density shows charge 
accumulation on the surfaces of the capacitor and local 
charge neutrality in the leads region as expected from a 
capacitor. 



B. Gold nanowires 

Metallic quantum point contacts (PCs) present con- 
ductance quantization at room temperature 5 - a property 
that has been predicted theoretically for many yearsi 5 -. 
Recently, Rodrigues et al. have shown that in a point 
contact the crystallographic orientation of the atomic 
tips forming the junction plays an important role in de- 
termining transport properties^. Therefore, a realistic 
theoretical description of the electronic transport in PCs 
must take into consideration the atomistic aspects of the 
problem. 

As an example we have performed calculations for a 
[100]-oriented gold quantum point contact (see inset of 
figure 0. A single gold atom is trapped at its equilib- 
rium position between two [100] fee pyramids. This is the 
expected configuration for such a specific crystal orien- 
tation, and the configuration likely to form in breaking 
junction experiments for small elongation of the junc- 
tion. This has been confirmed by atomic resolution TEM 
images^2iSi. In this case we have used LDA and a sin- 
gle zeta basis set for s, p and d orbitals. The unit cell 



of the extended molecule now contains 141 atoms (seven 
planes of the leads are included) and we consider peri- 
odic boundary conditions with 16 fc-points in the 2-D 
Brillouin zone. 




Figure 5: The transmission coefficient as a function of energy 
(upper panel) for a gold atomic point contact sandwiched be- 
tween two gold tips oriented along the [100] direction. In the 
lower panel the bandstructure for a monoatomic gold chain 
with lattice constant equal to the Au-Au separation in bulk 
gold. The inset shows a ball-and-stick representation of the 
atomic positions of the PC (the extended molecule) . 

In figure we present the zero-bias transmission coef- 
ficient as a function of energy. Recalling that the linear 
response conductance is simply G = 2e 2 /h T(Ep) (in 
this case we have complete spin-degeneracy) our calcula- 
tion shows one quantum conductance for this point con- 
tact. Interestingly the transmission coefficient is a rather 
smooth function T ~ 1 for a rather broad energy range 
around Ep. This means that the G = 2e 2 /h result is sta- 
ble against the fluctuations of the position of the Fermi 
level, which may be encountered experimentally. 

The large plateau at T ~ 1 indicates the presence of a 
single conductance channel for energies around and above 
E-p . This is expected from the bandstructure of a straight 
monoatomic gold chain with lattice parameter equal to 
the Au-Au separation in bulk gold (see figure (SJ)), which 
presents only one s band for such energy range. Therefore 
we conclude that the transport at the Fermi level is dom- 
inated by a single low-scattering s channel. Notably for 
energies 1 eV below Ep the transmission coefficient shows 
values exceeding one, which are due to contributions from 
d orbitals. In gold mono- atomic chains these are substan- 
tially closer to Ep than in bulk gold and participate to 
the transport. These results are in good agreement with 
previously reported calculations^^ and experimental 
d&taSLSM^. Additional examples of Smeagol's calcula- 
tions for PCs carried out by the authors can be found 
elsewhere in the literature 64,65 . 
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C. Molecular spin- valves 

The study of the I-V characteristics of magnetic sys- 
tems at the nanoscale is one of Smeagol's main goals. The 
most typical among spin devices is the magnetic spin- 
valve, which is obtained by sandwiching a non-magnetic 
spacer between two magnetic contacts. The direction of 
the magnetization in the two contacts can be arbitrarily 
changed by applying a magnetic held. The device then 
switches from a low resistance state, when the magneti- 
zation vectors in the leads are parallel to each other, to 
a high resistance state when the alignment of the mag- 
netizations is antiparallel. This is the giant magnetore- 
sistance e&eci£§£l, which is at the foundation of modern 
hard-disk reading technology. 

Traditional spin-valves use either metals or inorganic 
insulators as spacers. However a recent series of ex- 
periments have shown that organic molecules can serve 
the same purpose and a rather large GMR can be 
found§2i§2iZ2i£iiZ£. These experiments could lead to 
integrating the functionalities of molecules with spin- 
systems, and therefore have the potential to merge to- 
gether the fields of spin- and molecular-electronics. 

The calculation of the transport properties of molecu- 
lar spin- valves is a tough theoretical problem. It involves 
the computation of accurate electronic structures for 
magnetic surfaces, the charging properties of molecules, 
and the knowledge of the actual atomic positions. In 
a recent paper— we have demonstrated that molecules 
can efficiently be employed in spin-valves. Moreover 
we have shown that 7r-conjugated conducting molecules 
produce larger GMR then their insulating counterparts. 
Most of the effect is due to the orbital selectivity of 
the molecule/metal bonding, which in transition metals 
translates to a spin-selectivity. Here we further expand 
this concept and we demonstrate that the GMR can be 
tuned by molecular end-group engineering. 

The system under investigation is a 1,4-phenyl 
molecule attached to two fee Ni surfaces oriented along 
the (001) direction. The molecules are attached to the Ni 
hollow site through a thiol-like group where we use S, Se 
and Te as anchoring atoms. We consider collinear spin 
only and investigate the I-V characteristic assuming the 
magnetization vectors in the current /voltage contacts to 
be either parallel (P) or antiparallel (AP) to each other. 
The size of the GMR effect is expressed by the GMR ratio 
Rmr, which is defined as Rmr — (Ip — Iap)/Iap, with Ip 
(Jap) the current in the parallel (antiparallel) state. At 
zero bias, when all the currents vanish, we replace them 
with the conductances. 

We construct the unit cell of the extended molecule to 
include four Ni atomic planes on each side, for a total of 
forty Ni atoms. The basis set is critical and a single zcta 
for all the orbitals is not sufficient. Therefore we have 
used single zeta for H, C and S s orbitals, double zeta for 
Ni s, p and d, and double zeta polarized for C and S p 
orbitals. This basis gives us a Hamiltonian with over a 
thousand degrees of freedom. Finally the charge density 
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Figure 6: Transport properties for a 1,4-phenyl molecule at- 
tached to Ni (100) surfaces through a S group. The top 
panel shows the I-V characteristics for both the parallel and 
antiparallel alignment of the leads and the inset the corre- 
sponding GMR ratio. The lower panel is the transmission 
coefficient at zero bias as a function of energy. Because of 
spin-symmetry, in the antiparallel case we plot only the ma- 
jority spin. 



is obtained by integrating the Green's function over 50 
imaginary and 600 real energies. 

In figures H3 [7| and [S] we present the I-V character- 
istics, the zero bias transmission coefficient as a func- 
tion of energy, and the GMR ratio as a function of bias 
for the three anchoring situations (S, Se, Te). Clearly 
all the three cases show a large GMR, particularly for 
small biases. Interestingly the maximum GMR increases 
when going from S to Se to Te, and this is correlated 
with a general reduction of the total transmission and 
consequently of the current. Such a reduction is more 
pronounced in the case of antiparallel alignment of the 
leads and this gives rise to the increase in GMR. The 
origin of the drastic reduction of the transmission when 
changing the anchoring groups has to be found in the dif- 
ferent bonding structure. Since S, Se and Te all belong 
to the same row of the periodic table, the orbital nature 
of the bonding to the Ni surface is left unchanged, and 
so are the generic features of the transmission coefficient. 
However the bond distance goes from 1.28A to 1.48A to 
1.77A when going from S to Se to Te. This large increase 
in the bond distance is responsible for the reduction in 
transmission. 

The most relevant features of the transmission coeffi- 
cient can be understood in terms of tunneling through 
a single molecular stated. If we define t^{E) {t^{E)) as 
the majority (minority) spin hopping integral from one of 
the leads to the molecular state, then the total transmis- 
sion coefficient through the entire device in the parallel 
alignment will be simply T = T^T _|_ T ii = ( t T)2 + ^1)2 _ 
Here the total transmission coefficient for majority (mi- 
nority) spin is T TT = (t 1 ) 2 (T^ = (t 1 ) 2 ). Similarly in 
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Figure 7: Transport properties for a 1,4-phenyl molecule at- 
tached to Ni (100) surfaces through a Se group. The top 
panel shows the I-V characteristics for both the parallel and 
antiparallel alignment of the leads and the inset the corre- 
sponding GMR ratio. The lower panel is the transmission 
coefficient at zero bias as a function of energy. Because of 
spin-symmetry, in the antiparallel case we plot only the ma- 
jority spin. 



It is important to note that the large GMR ratio is 
ultimately due to the low transmission around E-p in the 
antiparallel case, which originates from the small trans- 
mission of the minority spins in the parallel case through 
the relation T n (E F ) oc ^T^ (E F )T±i(E F ). This is sur- 
prising since the density of states for minority spins at the 
Fermi level is rather large and one may expect substan- 
tial transmission. Moreover s-like electrons, which are 
weakly affected by the spin orientation, generally con- 
tribute heavily to the current regardless of the magnetic 
state of the device. Therefore, what does block the mi- 
nority spin electrons? 

A detailed analysis of the local density of states of the 
molecule attached to the leads^i reveals that the bonding 
of the molecule to the Ni surface is through Ni d and S p 
orbitals (the same is valid for Se and Te). The transport 
is therefore through hybrid Ni d-S p states, which in turns 
are spin-split due to the ferromagnetism of Ni. The cru- 
cial point here is that for the minority band these states 
end up above the Fermi level, and therefore do not con- 
tribute to the low bias transport. This is an important 
observation, since it demonstrates that orbital selectivity 
in magnetic systems can produce a spin-selectivity, and 
therefore magnetoresistance type of effects. 



the case of antiparallel alignment of the leads we have 
T = 2T n = 2T iT = 2tH l . In this simple description, 
which neglects both co-tunnelling and multiple scatter- 
ing from the contacts, T'^-(E) turns out to be a convolu- 
tion of the transmission coefficients for the parallel case 
T n oc VTTTTll. This type of behavior can be appreci- 
ated in figures H3 and |S] 
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Figure 8: Transport properties for a 1,4-phenyl molecule at- 
tached to Ni (100) surfaces through a Te group. The top 
panel shows the I-V characteristics for both the parallel and 
antiparallel alignment of the leads and the inset the corre- 
sponding GMR ratio. The lower panel is the transmission 
coefficient at zero bias as a function of energy. Because of 
spin-symmetry, in the antiparallel case we plot only the ma- 
jority spin. 



D. Nickel point contacts 

The transport properties of magnetic transition metal 
point contacts have been the subject of several recent 
investigations. Technologically these systems are attrac- 
tive since they can be used as building blocks for read 
heads in ultra-high density magnetic data storage de- 
vices. From a more fundamental point of view they offer 
the chance to investigate magnetotransport at the atomic 
level. Magnetic point contacts are effectively spin-valve- 
like devices, with the spacer now replaced by a narrow 
constriction where a sharp domain wall can nucleate*^. 
Therefore the magnetoresistance can be associated with 
domain-wall scattering and the MR ratio can be defined 
earlier. 

A simple argument based on the assumption that all 
the valence electrons can be transmitted with T ~ 1 gives 
an upper bound for the GMR of the order of a few per- 
cent (100% in the case of nickel). This however may be 
rather optimistic since one expects the d electrons to un- 
dergo quite some severe scattering. Indeed small values of 
GMR for Ni point contacts have been measured^. Sur- 
prisingly at the same time other groups have measured 
huge GMR for the same system^*2Si2L2i. Although me- 
chanical effects can be behind these large values^, the 
question on whether or not a large GMR of electronic 
origin can be found in point contacts remains. 

Therefore we investigate the zero bias conductance of a 
four atom long monoatomic Ni chain sandwiched between 
two Ni (001) surfaces (see figure 0I. This is an extreme 
situation rarely found in actual break junctions^?. How- 
ever an abrupt domain wall (one atomic spacing long) in 
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a monoatomic chain is the smallest domain wall possi- 
ble, and it is expected to show the larger GMR. For this 
reason our calculations represent an upper bound on the 
GMR obtainable in Ni only devices, and they also serve 
also as a test of the Smeagol capability for dealing with 
non-collinear spin. 




Figure 9: Schematic representation of the Ni point contact 
simulated. In the symmetric case the domain wall is located 
between the second and the third atom, while in the asym- 
metric it is placed between the third and the fourth. The 
direction of the current is from 1 to 4 for positive bias. 

In this calculation we use a double zeta basis set for 
s, p and d orbitals and consider finite leads (no peri- 
odic boundary conditions are applied) with either four 
or five atoms in the cross section. We then investigate 
two possible situations. In the first one we place the do- 
main wall symmetrically with respect to the leads, i.e. 
between the second and the third atom of the chain. In 
the second (asymmetric) the domain wall is positioned 
between the third and the fourth atom. Furthermore 
we perform spin-collinear and spin-non-collinear calcula- 
tions for both cases. Interestingly all our non-collinear 
calculations always converge to a final collinear solution. 
This confirms expectations based on simple s-d models, 
suggesting that the strong exchange coupling between 
the conduction electrons and those responsible for the 
ferromagnetism, stabilize the collinear state if the mag- 
netization vectors of the leads are collinear. 

In figure^JIwe present the transmission coefficient as a 
function of the energy for both the symmetric and asym- 
metric case and the parallel state. For collinear calcula- 
tions the contributions from majority and minority spins 
are plotted separately, while we have only one transmis- 
sion coefficient in the non-collinear case. Clearly in all 
cases the non-collinear solution agrees closely with the 
collinear one, i.e. _ T c T ollinoai . + = T non _ conincar . 

This is expected since the final magnetic arrangement of 
the non-collinear calculation is actually collinear, and it 
is a good test for our computational scheme. 

Turning our attention to the features of the transmis- 
sion coefficient it is evident that at the Fermi level T in 
the parallel state is larger than that in the antiparallel. 
This difference however is not large and the GMR ratio 
is about 60% with little difference between the symmet- 
ric and asymmetric domain wall. This is mainly due to 
the much higher transmission of the un-polarized s elec- 
trons compared with that of the d. Note that the conduc- 
tance approaches 2e 2 /h for energies approximately 0.5 eV 
above the Fermi level. For such energies in fact no d elec- 
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Figure 10: Transmission coefficient as a function of energy 
for the nickel quantum point contacts of figure The right- 
hand side panels (a) are for collinear calculations and the 
left-hand-side (b) are for non-collinear; (1) parallel state, (2) 
antiparallel with symmetric domain wall, (3) antiparallel with 
asymmetric domain wall. Note that in the non-collinear case 
we do not distinguish between majority and minority spins. 
In panel (a2) majority and minority spins are degenerate. 



trons contribute to the density of states of both the spin 
sub-bands, and only s electrons are left. These are then 
transmitted with T ~ 1 as in the case of Au chains in- 
vestigated previously. 

The crucial point is that the contribution of the s 
electrons is also large at the Fermi level. This results 
in a poorly spin-polarized current at low bias and con- 
sequently in a small GMR, in agreement with other 
calculations^ 1 *^. In conclusion our finding rules against 
the possibility of large GMR from electronic origin in Ni 
point contacts. However the presence of non-magnetic 
contamination (for example oxygen) may change this pic- 
ture radically. 



E. H2 molecules joining platinum electrodes 

The aim of this section is to show how strongly the 
transmission coefficients may depend on the leads cross 
section, whenever d electrons are close to Ep. As an ex- 
ample, we present results for the H2 molecule sandwiched 
between fee Pt(001) leads, and compare leads of differ- 
ent sections with extended leads. These are obtained by 
applying periodic boundary conditions and a sampling 
over the fc-points along the direction orthogonal to that 
of transport. 

The conductance of an H 2 molecule sandwiched 
between platinum electrodes has been extensively 
studied 64 i 83 i 84 i 85 i 86 . Experimentally 8 ^ it has been found 
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that the inclusion of hydrogen gas into the vacuum cham- 
ber produces a dramatic change in the conductance his- 
tograms of platinum, which change from a structure with 
a broad peak at 1.5 Go to a structure with a sharp peak 
at 1 Go- This resonance has been attributed to the con- 
ductance through a single molecule, which bridges both 
leads lying parallel to the current flow. This explanation 
has been confirmed by theoretical calculations^^iSi. 
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Figure 11: Transmission coefficients for an H2 molecule sand- 
wiched between fee platinum leads near the equilibrium dis- 
tances (~ 9.5 and 11 A). The left hand side (a) corresponds 
to the configuration where the molecule lies parallel to the 
current flow and the right hand side (b) to the configuration 
where the molecule lies perpendicular. The leads are made of 
alternating slabs of 4-5 atoms (1) and 9-12 atoms (2) without 
periodic boundary conditions along the perpendicular direc- 
tions (xy), and 9-9 atoms with periodic boundary conditions 
along xy (3). In the last case the dashed and continuous lines 
have been obtained with 4 and 12 k points, respectively. 

In our calculations the H2 molecules is located either 
parallel or perpendicular to the current flow. A detailed 
description of the geometric configuration and the re- 
sults can also be found in reference |6J| . We use a double 
zeta polarized basis set for platinum s, p and d orbital, 
a double zeta for the hydrogen s electrons and the LDA 
functional. As a first step we employ finite cross sec- 
tion leads along the transversal directions, composed of 
alternated planes containing 4 and 5 atoms each. The 
resulting transmission coefficients show many peaks and 
gaps throughout all the energy range and particularly 
sharp variations around the Fermi energy, as can be seen 
in figureslTDfal) and (bl). When thicker slabs composed 
of alternating planes of 9 and 12 atoms are employed the 
results do not improve and the large oscillations still re- 
main, as shown in figures im fa2) and (b2). It is apparent 
from these figures that while T(E) shows a long plateau 
at positive energies, it presents strong oscillations at the 



Fermi energy and, therefore, it is uncertain to infer the 
conductance of the junction from T(Ep). 

This is in stark contrast with the case of gold, where 
the d- levels lie below Ep, and T(E) is smooth regardless 
of the size of the leads cross section. For platinum the 
presence of d-states at the Fermi energy opens minibands 
and minigaps, which translate into strong oscillations in 
T(E ~ Ep). These minibands and minigaps arise from 
interference effects of the d-states along the transverse 
direction. Consequently, oscillations in T(E) should dis- 
appear when bulk electrodes are used. Indeed, this is 
what we find when slabs made of 3 x 3 atomic planes and 
periodic boundary conditions are employed, as shown in 
figures ITTI (a3) and (b3). We moreover show how T(E) 
converges when the number of transverse k points is in- 
creased from 4 to 12. Although some small variations and 
peaks still remain when 4 k points are used, the trans- 
mission at the Fermi level is essentially converged. Note 
that the parallel case has T ~ 1 for a long range of ener- 
gies around Ep, which remains essentially unperturbed 
for small variations of the coordinates or the distance 
between the electrodes. This explains the sharp peak 
observed in the experimental conductance histograms 83 . 

In view of the above calculations we can therefore con- 
clude that the use of bulk electrodes, characterized by 
periodic boundary conditions along the perpendicular di- 
rections and k points, is mandatory in order to avoid 
oscillations in the transmission coefficients. Otherwise 
the presence of strong variations and minigaps can give 
unphysical solutions for systems with open d shells. 



IV. CONCLUSIONS 

We have presented a description of our newly devel- 
oped non-equilibrium Green's function code Smeagol. In 
the present version Smeagol uses the DFT implemen- 
tation contained in Siesta as the underlying electronic 
structure method. However the code has been developed 
in a modular and general form and can be easily com- 
bined with any electronic structure scheme based on lo- 
calized orbital basis set. The core of Smeagol is our new 
algorithm for calculating the surface Green's functions 
of the leads, which combines generalized singular value 
decomposition with decimation. This results in an un- 
precedent numerical stability for a quantum transport 
code, and in the possibility of drastically reducing the 
number of degrees of freedom in the leads. In this way 
large current /voltage probes with complicated electronic 
structure can be tackled. 

We have also presented a selection of results obtained 
with Smeagol. These range from simple tests for the elec- 
trostatics to an analysis of the GMR in molecular spin- 
valves, and demonstrate the SmeagoFs ability to tackle 
very different problems. 
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Appendix A: STABLE ALGORITHM FOR THE 
EVALUATION OF THE SELF-ENERGIES 

1. The "A'i problem" 



The method presented in section Hi Dl to calculate the 
leads Green's functions depends crucially on the fact that 
the coupling matrix between principal layers K\ = Hi — 
ESi is invertible and not ill-defined. However this is not 
necessarily the case since singularities can be present in 
Ki as the result of poor coupling between PLs or because 
of symmetry reasons. Note also that since Ki = Hi—ESi 
the rank of K\ may also depend on the energy E. 

We now give a few examples illustrating how these sin- 
gularities arise. Let us consider for the sake of simplic- 
ity an orthogonal nearest neighbour tight-binding model 
with only one s-like basis function per atom. In this case 
K\ = Hi is independent from the energy. In figure IT21 
we present four possible cases for which Hi is singular. 
In the picture the dots represent the atomic position, the 



(a) 





Figure 12: Four different structures for which Hi is singular: 
(a) lack of bonding, (b) supercell, (c) over-bonding, (d) odd 
bonding. Each black dot represents an atom and each line a 
bond. The dashed boxes enclose a principal layer. 

lines the bonds and the dashed boxes enclose a PL. All 



the bonds are assumed to have the same strength, thus 
all hopping integrals 7 are identical. 

In the first case (figure ^Ji) the PL coincides with the 
primitive unit cell of the system and therefore it is the 
smaller principal layer that can be constructed. However 
since every second atom in the cell does not couple with 
its mirror in the two adjacent cells Hi has the form 



Hi 



7 




(Al) 



and therefore is singular. This is the case of "lack of 
bonding" between principal layers. It is the most com- 
mon case and almost always present when dealing with 
transition metals, since localized d shells coexist with de- 
localized s orbitals. 

Figurelulb presents a different possibility. Here the PL 
is a supercell constructed from two unit cells and every 
atom in the PL couples with atoms located in only one 
of the two adjacent PLs. In this case 



Hi = 





7 



(A2) 



which is again singular. Clearly in this specific case one 
can reduce the principal layer to be the primitive unit cell 
solving the problem (Hi become a scalar 7). However 
in a multi-orbital scheme the supercell drawn may be 
the smallest PL possible and the problem will appear. 
Again this is a rather typical situation when dealing with 
transition metals. 

The case of "overbonding" is shown in figure 112b . 
Again the PL coincides with the primitive unit cell, but 
now every atom in the PL is coupled to all the atoms in 
the two adjacent PLs. In this case 



Hi 



1 1 
1 1 



(A3) 



which is not invertible. 

Finally the "odd bonding" case is presented in figure 
I12H . Also in this case the PL coincides with the primitive 
unit cell, however the upper atom in the cell is coupled 
only to atoms in the right nearest neighbour principal 
layer. The Hi matrix is then (we label as "1" the upper 
atom in the cell) 



Hi = 



7 
7 



(A4) 



i.e. it is singular. Clearly the above categorization is 
basis dependent, since one can always find a unitary ro- 
tation transforming a generic Hi in a new matrix of the 
form of equation (|A4|) . 



2. Finding the singularities of K\ 

We now present the first step of a scheme for reg- 
ularizing Ki, and indeed the whole Hamiltonian and 
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overlap matrix, by removing their singularities. In the 
cases of "lack of bonding" , "supercell" and "odd bond- 
ing" presented in the previous section the singularities of 
K\ = Hi were well defined since an entire column was 
zero. However more generally, and in particular in the 
case of multiple zetas basis set, K\ is singular without 
having such a simple structure (for instance as in the 
"overbonding" case). This is the most typical situation 
and a method for identifying the singularities is needed. 

The ultimate goal is to perform a unitary transforma- 
tion of both Tt and S in such a way that the off-diagonal 
blocks of the leads Hamiltonian and overlap matrix (Hi 
and Si) assume the form 



N-R R 

n [0 , A] 





• 


• A hN 


-R+l 





• 


■ A 2 ,n 


-R+l 





• 


■ A 3: n 


-R+l 



■■ A 1<N \ 

■ ■ A 2 ,n 

■ ■ A 3 M 



\0 



At 



^n.n-r+i ■ ■ ■ An n / 

(A5) 

i.e. they are N x N block matrices of rank R, whose 
first N — R columns vanish. In this form the problem is 
re-conducted to the problem of "odd bonding" presented 
in the previous section. 

This can be achieved by performing a generalized sin- 
gular value decomposition (GSVD) 87 . The idea is that a 
pair of N x N matrices, in this case Hi and Si, can be 
written in the following form 



Hi = UAi [0,W}Q\ 
Si = VA 2 [0,W] Q\ 



(A6) 
(A7) 



with U, V and Q unitary N x N matrices and W being a 
Rx R non-singular triangular matrix where R is the rank 
Hi ' 



of the 2N x N matrix 



Si 



(R < N). The matrices A : 



and Ao are defined as follows: 



Ai 



K L 




(A8) 



A 2 



N- 



K L 

o a 
o o 



(A9) 



where L is the rank of Si, K + L = R, Ik is the K x K 
unit matrix and C and C are matrices to determine. 

Clearly the two matrices Hi and Si have the two com- 
mon generators W and Q. Then, one can perform a 
unitary transformation of both Hi and 5i by using Q, 
obtaining 



H^ = Q^HiQ = Q^UKi [0,W] 
= Q f SiQ = Q f VA 2 [0,W]-- 



N-R R 

n [ , Hi ] 

N-R R. 



N 



[0 , Si] 



(A10) 
(All) 



Here Hi and Si are the N x R non-vanishing blocks of 
the GSVD transformed matrices Hi and Si respectively. 

In an analogous way the same transformation for H\, 
S\, H and So leads to 



H 



q?h\q 



sf = q^s\q 



N-R 
R 



N-R 
R 



N 



N 

' 



= Q^H Q , 
S Q = Q+SoQ , 



(A12) 

(A13) 

(A14) 
(A15) 

are not nec- 



where the transformed matrices H® and Sq 
cssarily in the form of equation l|A5(l . 

We are now in the position of writing the final uni- 
tary transformations for the total (infinite) Hamiltonian 
H and overlap S matrices describing the whole system 
(leads plus extended molecule). These are given by 
Q^HQ and Q^SQ with the infinite matrix Q defined as 



Q 



V 



\ 



Q . 
Q 
■ I M 
. . 



. 

Q 
Q 



(A16) 



/ 



where Im is the M x M unit matrix. Note that this 
unitary transformation rotates all the Hi matrices (the 
Si matrices in the case of S), but it leaves Hm (Sm) 
unchanged. Finally the matrices coupling the extended 
molecule to the leads transform as following 



n LM 

H Q 

H Q 

RM 



H 



MR. 



HmlQ , 
HmrQ , 



(A17) 



and so do the corresponding matrices of S. 



3. Solution of the "A 1 problem" 

Now that both H. and S have been written in a conve- 
nient form we can efficiently renormalize them out. The 
key observation is that the two (infinite) blocks describ- 
ing the leads have now the following structure (the iS 
matrix has an analogous structure and it is not shown 
here explicitly) 
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H Q - 

L/R — 



Q^HqQ Q^H.Q 
Q^H-xQ Qt^ g Qt^g 

o gt^-ig q^oO 



Q^HoQ [0 Hi] 












5 


I) 





[0 ff! 

















(A18) 



where the matrices D, B and C are respectively R x R, 
N x (N -R) and (N - R) X (N - R). 



Note that the degrees of freedom (orbitals) contained 
in the block C of the matrix Hq — Q^HqQ couple to 
those of only one of the two adjacent PLs. This situation 
is the generalization to a multi-orbital non-orthogonal 
tight-binding model of the "odd bonding" case discussed 
at the beginning of this appendix (figure \12\ i) . These 
degrees of freedom are somehow redundant and they will 
be eliminated. We therefore proceed with performing 
Gaussian elimination^, (also known as "decimation") of 
all the degrees of freedom associated to all the blocks C. 

The idea is that the Schrodinger equation Q'[H — 
ES]Q^! = can be re-arranged in such a way that a 
subset of degrees of freedom (in this case those associ- 
ated to orbitals in a PL that couple only to one adjacent 
PL) do not appear explicitly. The procedure is recur- 
sive. Let us suppose we wish to eliminate the l-ih row 
and column of the matrix KP = Q)[H — ES]Q. This can 
be done by re-arranging the remaining matrix elements 
according to 



K9 a K9 



ICQ, 



(A19) 



The dimension of the resulting new matrix ("1" 
indicates that one decimation has been performed) is re- 
duced by one with respect to the orig inal/C Q . This pro- 
cedure is then repeated and after r decimations we obtain 
a matrix 



= k, q 



0) 



K9 



(r-1) 



(A20) 



Let us now decimate all the matrix elements contained 
in all the sub-matrices C. We obtain a new tridiagonal 
matrix "oo" means that an infinite number of 



decimations have been performed) of the form 



K, Q 



(oo) 



et 
o 



V 



A 

et 
o 



e 

A 

ti 







Ti 
Di 
K 



Q 

ML 







K Q 

n LM 

K M 

©RM 







©MR 

D 2 

et 
o 



o 
e 

A 

et 



o 
e 

A 



H 



LM 



ES 



LM' 



K- 



Q 



ML 



h: 



Q 



ML 



(A21) 
ES ML and 



where K£ u 

Km = Hm — ESm- The crucial point is that the new ma- 
trix K.Q(°°) is still in the desired tridiagonal form, but now 
the coupling matrices between principal layers © are not 
singular. These are now RxR matrices obtained from the 
decimation of the non-coupled degrees of freedom of the 
matrices (the C blocks). Moreover the elimination of 
degrees of freedom achieved with the decimation scheme 
is carried out only in the leads. The electronic structure 
of these is not updated during the self-consistent proce- 
dure for evaluating the Green's function, and therefore 
the information regarding the decimated degrees of free- 
dom are not necessary. In contrast the degrees of freedom 
of the scattering region are not affected by the decimation 
or the rotation. Therefore the matrix Km is unaffected 
by the decimation. 

In the decimated matrix JC®^ 00 ^ new terms appear (D\, 
D 2 , T\ and ©mr)- These arise from the specific structure 
of the starting matrix [H — ES] Q and from the fact 
that the complete system (leads plus scattering region) is 
not periodic. In fact assuming that j is the last principal 
layer of the left-hand side lead and I is the first layer of 
right-hand side lead, the decimation is carried out up to 
j — 1 to the left and starts from I to the right of the scat- 
tering region. This allows us to preserve the tridiagonal 
form of K. and at the same time to leave Km unchanged. 
A schematic picture of the decimation strategy is illus- 
trated in figure l*Hfl 

In practical terms all the blocks of the infinite matrix of 
equation (|A21fl can be calculated by decimating auxiliary 
finite matrices. In particular 

1. A, © and D 2 are calculated by decimating both the 
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(a) 



(b) 



(c) 



Figure 13: Schematic representation of the decimation strat- 
egy for the rotated /C matrix . Every symbol (dots, boxes 
..) represents a collection of degrees of freedom (a matrix 
block) and every line the coupling, (a) Original structure af- 
ter the rotation Q. In the periodic leads the upper black dots 
represent the blocks C of the matrix of equation 1A18L The 
large white rectangular box represents the scattering region, 
(b) The degrees of freedom marked with the red crosses are 
decimated, (c) Final structure after decimation. The new 
white symbols represent the leads degrees of freedom of the 
principal layers adjacent to the scattering region as they ap- 
pear after the decimation. 



C matrices of the following finite 2N x 2N matrix 
\ 



C B 
St D 



[OK,] 












V 


M. 


(S 


2) 



d 2 e 

0t A 



(A22) 



where K\ = H x — ESi and 



C B 
St D 



(A23) 



2. Di and T\ are calculated by decimating only the 
upper C matrix of the same finite 2N x 2N matrix 



( ( C B 



[0 Kt 












V 




(S 


2) 



D 2 Ti 



(A24) 



where D x is N x N, while T x is R x N. 

3. 9mr is a M x R matrix obtained by decimating 
the C block of the following (N + M) x (N + M) 
matrix 







M 

Q 

RM 




0m ©MR 



(A25) 



where 0m is the M-dimensional null matrix. 

Finally we are now in the position of calculating 
the self-energies. These are obtained from the surface 



Green's functions for the rotated and decimated leads 
(specified by the matrices A and 0) and have the follow- 
ing form 



Q 

LM 



and 



Si? = ©mr [Gr 1 - (D 2 - A)] 1 e 



RM 



(A26) 



(A27) 



Clearly our procedure not only regularizes the algorithm 
for calculating the self-energies, giving it the necessary 
numerical stability, but also drastically reduces the de- 
grees of freedom (orbitals) needed for solving the trans- 
port problem. These go from (the dimension of the 
original Hi matrix) to R (the rank of Hi). Usually 
R -C N and considerable computational overheads are 
saved. 

Finally it is important to note that usually the rank R 



of 



Hi 
Si 



is not necessarily the same of that of 



St 



(i?'). If R' < R the GSVD transformation must be per- 
formed over the matrices H[ and S[. The procedure is 
similar to what described before but the final structure 
of the matrix KP is somehow different, and so should be 
the decimation scheme. 



Appendix B: THEORETICAL DESCRIPTION 

1. Brief reminder of Keldish formalism 

The electronic part of the system we consider is de- 
scribed in a general way by the following Hamiltonian 88 

H{r u r 2 ,t) = H {ri,t)5{fi - r 2 ) + H x {r u r 2 ) (Bl) 

where H x accounts for the Coulomb interaction among 
electrons and H stands for all one-particle pieces of the 
Hamiltonian, 



H (f,t) = ff k in(f) + H ci {r) + V cxt (r, t) 



(B2) 



i/kin and H el are respectively the kinetic energy and the 
Coulomb interaction between electrons and nuclei, and 
V ex t is an external electrostatic potentials applied to the 
system at t = 0. 

The electronic system evolves according to the time- 
independent Hamiltonian H(to = 0~) for all negative 
times, and uses to as the synchronization time for all pic- 
tures. Moreover we assume that before V cx t is switched 
on the system of interacting electrons is in thermody- 
namic equilibrium at a chemical potential /^o- We then 
prepare the density matrix p at to also. Therefore, ex- 
pectation values of observables 

(6(t)) = Tr{p s (t)6 s (t)} = Tr{p H (t Q )d H (t)} 

= _L Tr{e -P { H{ to )-„ N)QH^ (B3) 
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are described in terms of the density matrix of an inter- 
acting electron system at equilibrium. 

These expectation values may be evaluated by using 
perturbation theory Wick's theorem may be applied 
only to ensembles of non-interacting electrons. In order 
to take advantage of it, we must use a non- interacting 
density matrix such as 



Po(io) 



-0{H o (t o )-^ o N) 



(B4) 



We then define the Hamiltonian in the Schrodinger pic- 
ture in the Keldish contour of Fig. ^] of the complex r 
plane as follows 



K s (r) = 



h s (t) t e c H 

H s (0)-ij, q N t e c v 



(B5) 



and analogous expressions for its noninteracting and in- 
teracting pieces, K§(t) and Kf{r). We notice that the 
time variable is doubled valued along the real axis. We 
must therefore distinguish whether any real time lies on 
the upper (t + ) or lower (t~) branch. 



Rex 



C, 



ImT 



Figure 14: Keldish contour c in the complex r plane. The 
imaginary time path is called ci or cv in the text. The seg- 
ment lying below (above) the real-time axis is called ci icz). 
The time loop C2+C3 is called ch- 

We then define Heisenberg (H) and Dirac (D) pictures 
in the Keldish contour whereby evolution of operators is 
provided by 



which are evaluated in the non-interacting ensemble de- 
scribed by pq. 

The Green function of the system, 



G(1,2) = <T C ^(1)^ T (2)> 



(B7) 



is a tool to compute the physical response of the system. 
Thus the electron charge and current densities are simply 



(n(l)> 



-ieG(l+ 1" 
eh 



(B8) 



Oil)) = -7T- V 1 -V 2 + 2zA(l) G(l+2)| 2=1 - 

and depend on the applied the external potential. Here, 
ip , %jr denote creation and annihilation operators, and 
(i) = (fi,U). 

The Green's function can be proven to satisfy the fol- 
lowing Dyson equation 

[id tl -H (l)] G(l,2) = «5(1,2) + (E®G)(1,2) 
[-id t2 -H {2)\ G(l,2) = *(1,2) + (G®E)(1,2) 

where all time variables run through the Keldish contour 
c and we follow the conventional shorthand 



(4® B)(l,2) 



dx 3 dt 3 A(l,3)B(3,2) 



(B9) 



Connection to TDDFT 



We now define a fictitious system of non-interacting 
electrons described by the Hamiltonian^ 



H s (f,t)=H Mn (r) + V s (r,t) 



(BIO) 



where we assume that H s is constant and that the sys- 
tem is in thermodynamic equilibrium at chemical poten- 
tial /j, s , all for negative times. We also assume that H s 
may be exactly diagonalized. We do not split a perturb- 
ing piece from the Hamiltonian and therefore, Heisenberg 
and Dirac pictures coincide. We again take to as the syn- 
chronization and density-matrix preparation time. The 
density matrix 



Ps(to) 



-0{H s {t o )-,x 3 N) 



(Bll) 



and the time-evolution operator on the Keldish contour 



o H (t) 


— T c e 




dr' 


KS{T,) 6 s {t) 


o D (t) 


= T c e 




dr' 


K ° {T,) d s (t) 




= T c e 


~ l Io+ 


dr' 





U S (0 + ,Q-)=T C e 



(B12) 



'6*(t) =T c U{0 + ,0-)O s (t) 

, _ ~ s completely determine the expectation value of observ 
) ,0 )0 (t) a jjj eg f ^jjjg fictitious system 

-i f o ° + *SK?{r>) 6 H {t) = Tc Vl{ Q+ :{r) d*(t) { o(t)) s = Tr{p s s (t) 6 s (t)} = Tr{pf (< ) 6 H (i)} 

where T c is the time ordering operator on c. Expecta- 
tion values of operators are then given by the statistical 
averages 

Zo 



The Green's function G s (l,l'), whose time variables 
also lie on c, satisfies the equation of motion 



(o(t)) 



z 



Tr{ Po V l (-if3^ + )0{t)}. (B6) 



[id tl 
-idu 



H s (l)} G s (l,2) 
H s (2)} G s (l,2) 



8(1,2) 
8(1,2) 



(B13) 
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and provides the physical response of the system 



n s [V s (l)} = -*eG a (l+ I") 



2 m 



V 1 -V 2 + 2zA(l)j G s (l+,2)| 2=1 
We define the action functionals 

R[V(x)] = UnTr \u(-i0,O) 
R s [V s (x)] = HnTr[{7 s H/3,0)j (B14) 

whose functional derivative with respect to the external 
potential provides an alternative way to calculate the 
electronic density. We adjust the potential V s so that 
the densities of the two systems be equal, 



SR[V(x)} SR[V a (x)] 



SV(x) 



5V s {x) 



= (n(x)) 



(B15) 



Notice that we do not require that the current densities 
of the two systems be equal. 

We now perform a Legendre transform to find the new 
actions 



S[n] 
S s [n] 



-R[V[n}\ + J dxn(x)V(x) 
-R s [V s [n}] + / dxn(x)V s (x) (B16) 



such that 



SS[n(x)} 

Sn(x) 
SS s [n(x)} 
Sn(x) 



= V(x) 
= V s (x) 



(B17) 



We define an exchange-correlation functional 

S s [n] = S[n] + \Jdx [Vn(x) n(x) + S xc [n}} (B18) 

whose functional derivative with respect to the density 
provides with a key relationship between the external 
potential of the actual and fictitious systems, 



V s (x) = V(x) + Vu(x)+V xc (x) 

ss xc 



V xc [n(x)] 



Sn(x) 



(B19) 



Comparing the equations of motion of G and G s , we 
find an explicit relationship between the Green's func- 
tions 

G = G s + G s ® ( £ - Vk - ) ® G (B20) 
and the Sham-Schluter equation for the Self-energy 

G a ®(V K + V xc )®G = G a ®Tl[C\®G (B21) 



Iterating these equations once means equating both 
Green's functions, G = G s , which in turn implies ap- 
proximating )byj s . The resulting integral equation for 
£ has 

X[G](1,2) = (Vk[n(l)] + KcKl)])<5(l,2) (B22) 

as a trivial solution. This simple approximation has the 
virtue that charge is conserved. This can be seen in two 
alternative ways. First, the self-energy can be written as 
a <£>[G] -derivable function. Second, the fictitious system 
satisfies the continuity equation by construction. Subse- 
quent iterations improve the physical content of G, but 
we have used in our code this lowest-order approxima- 
tion. 

The Green's function at two different times G(t, t') can 
be viewed as the matrix element (t\G\t'), sandwiched be- 
tween two time states of the whole set {\t}} of times in 
the Keldish contour. We split the contour into the three 
pieces c\, C2, C3 of Fig. 2, and define the corresponding 
three time subsets. The Green's function can then be 
represented in matrix form as 



G s (t,t') 



G n (M') G 12 (t,t') G 13 (t,t') 

G 21 (t,t r ) G 22 (t,t r ) G 23 (t,t r ) 

G 31 (t,t') G Z2 {t,t r ) G 33 (t,t') 

G c (t,t') G<(t,t') G 13 (t,t') 

G > {t 1 t') G ac (t,t r ) G 13 (t,t r ) 

G 31 (t,t r ) G 31 (t,t r ) G 33 (t,t r ) 



The physical response of the system is therefore encap- 
sulated in the lesser Green function, 



n 



,(i)] = -*G<(i + ,n 

n 



= ~2^ - V 2 + 2i^(l)J G<(1+,2)| 2=1 - 

There are only five independent G y out of the seven 
matrix elements displayed above. A partial reduction to 
six matrix elements in the Green's function is achieved 
by the non-unitary transformation 

G R G K V2G 13 
G s (t,t') = LT 3 G s (t,t')L^ = \ G A 

V2G 31 G 3 3 

Both matrices satisfy the matrix equations of motion 

[id tl -H,(l)] G s (l,2) - <5(ri - r 2 ) 5(n - r 2 ) 
\id tl -H S {1)} 6,(1,2) = 5{f 1 -f 2 )S{n-T2) 

where the time delta functions are defined as 

S(h - 1 2 ) 

5,5= \ ±5(h-t 2 ) I (B23) 
5{n - t 2 ) 

with the menos (plus) sign corresponding to the hat 
(check) delta function. We shall drop the "s" subindcx 
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from now on, since all the discussion that follows is valid 
only for the ficticious system of non-interacting electrons. 

Physical response is customarily written in terms of 
G R A and G < . Since there is no linear transformation 
which allow to group them in one common matrix, one 
uses Langreth's rules to find relationships among them. 

3. Equations of motion in the localized wave 
function basis 

The eigenstates of the system can be obtained by 
expanding them in the basis of non-orthogonal states, 
*(n,r,t) = Cin(n,t) ip^f - Ri), in terms of which 
the Schrodinger equation reads 

[iS^^dt.-H^j^h)} Cj V (n,h) = (B24) 

Alternatively the equation of motion for either check 
or hat Green functions are 

[iSip k\ d tl - Hin kx(ti) } Gk\ 3 u(ti,t 2 ) = Si 3 6nv5(ti—t2) 

(B25) 

It is advantageous to perform a change of time vari- 
ables from ti, £ 2 to T = l/2(ti + t 2 ),t = ti - t 2 . Then 
the Green's functions can be written as 

Giu ju(ti,t2) = Gi„ jv(T,t) = ( — — Gin j V (T, E) e lEt 

J 27T 

(B26) 

The electron charge and current densities are found 
from 

n(r,T) = n itl jv {f) p itl jv {T) 

%H jv 

W, T) = J2 3* J- ) i» ( T ) (B27) 

in jv 

where we have introduced the density matrix 

Pinjv{T) = J ^L G < JU (T,E) (B28) 

and 

rii/i j v {fi,f 2 ) = e ^(ri - Ri) ip v {r 2 - Rj) 

i e H / -* -* -* \ 

JipjAr) = -7) — ( V n - v ?2 +2iA(f 1 )J n VJU (r u r 2 )\ 

The electric current through a given surface S is ob- 
tained by integrating the current density over such sur- 
face 

I = dS-y(r,h)= ^2 PinjAti) / dS-j ifijv 



— Pi'n' j'v Hypi ji v i (B29) 

i'H 1 j'v' 

where only those bonds (i'/j.j'u') pierced by the surface 
contribute to the summation. 

If the system is in thermodynamic equilibrium, the 
population of electrons does not depend on time, and fol- 
lows the Fermi-Dirac distribution function /(e). Thereby 
the lesser Green function can be written in terms of the 
retarded one as 

Gf^jviE) = f(E) [G^ jv {E) - Gf^ jv {E)] = 
= -2if(E)lm [G% jv {Ej\ = 
= 2nif(E)^2c ill (n)c%(n)6(E-e n ) 

n 

where e n are the eigenvalues of the Hamiltonian. There- 
fore, the density matrix can be obtained from the wave 
function coefficients by just diagonalising the Hamilto- 
nian, 

/W = ~f dEf(E)Im [G? ,•„(£)] = 

= ]T c v (n)4»/(e n ) (B30) 

n 

4. The extended molecule setup 

We now wish to partition the Green functions accord- 
ing to the system setup of Fig. 1, where the left and right 
leads remain in thermodynamic equilibrium defined by 
Ml/r = E F ± eV/2 at all times. The extended molecule 
is also in thermodynamic equilibrium for times t < 0. 
The Hamiltonian for negative times 

/ Hl + eV/2 «Sl \ 

H(t < 0) = h{t) = H M 

\ H R - eV/2 S R J 

(B31) 

serves to define the reference equilibrium hat and check 
Green functions, that satisfy the equations of motion 

[iSi^kxdti - hi^kxih)] Gkx,jv{ti,t 2 ) = Sij S^^ S(ti-t 2 ) 

(B32) 

For instance, the equation of motion for the retarded 
Green function in frequency domain is just 



21 

















r<oR 




1 = 


i 




s 

















(e+ - eV/2)S L - H L \ / S£ R 

e+S M -H M ] I G$ | = | I u | , (B33) 

(e+ + eV/2)S RU - H RM 

while the lesser Green function is 

(Gl A (E)-Gl R (E)) HE-eV/2) 
g°<(E)=\ G°<(E) ) (B34) 

(g K \E)-g^{E)) f{E + eV/2) 



The extended molecule is contacted by the electrodes 
at time t — 0, through the Hamiltonian matrix elements 

/ H LM \ 

V ext = Hml H MR F(t), (B35) 
V o H RM o / 

where F(t) is zero for negative times, and 1 for times 
larger than a certain characteristic time tm ■ The pertur- 
bation V ex t drives the core of the extended molecule out 
of equilibrium for positive times, by populating it with 
a distribution of electrons that does not follow Fermi- 
Dirac statistics. The distribution function of the non- 



contacted molecule gl^(E) is completely washed out, but 
the density-matrix of the system can still be determined 
from the equations of motion of the lesser and retarded 
Green functions. 

We seek to solve the equations of motion for times 
t ^> tm where all transient effects have vanished. This 
means, first, that G 1 ' 3 = G 3 ' 1 = 0, so that the matrix hat 
and check Green functions are block-diagonal; second, 
that the Hamiltonian is simply 7i; and third, that Green 
functions do not depend on the time variable T. 

The retarded Green function is, simply, 



(e+ - eV/2) S h - H h e+S LM - H LM 

e + 5 RM - Hrm (e H 





e + <?MR — ^MR 

+ eV/2)S R -H R 



G R G R 



r?R r R c R 



RL y RM 



LR 
R 

MR 
R 
R 



g n 














• 










, (B36) 



The lesser Green function is better expressed in terms 
of the retarded and advanced, and the reference equi- 
librium lesser Green function, by using Langreth's rules, 
as 

g< = g R (g 0R )-i g°< (g^y 1 g A (B37) 



Straightforward matrix algebra then leads to Eqs. I jlfijl 
and CD). 
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